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The mathematical structure of Sudoku puzzles is akin to hard constraint satisfaction problems lying at the 
basis of many applications, including protein folding and the ground-state problem of glassy spin systems. 
Via an exact mapping of Sudoku into a deterministic, continuous-time dynamical system, here we show that 
the difficulty of Sudoku translates into transient chaotic behavior exhibited by this system. We also show 
that the escape rate k, an invariant of transient chaos, provides a scalar measure of the puzzle's hardness that 
correlates well with human difficulty ratings. Accordingly, f/ — — logio a- can be used to define a 
"Richter"-type scale for puzzle hardness, with easy puzzles having 0 < y/ < 1, medium ones 1 < ij^2, hard 
with 2 < // < 3 and ultra-hard with > 3. To our best knowledge, there are no known puzzles with ii> \. 



In Sudoku, considered as one of the world's most popular puzzles', we have to fill in the cells of a 9 X 9 grid with 
integers 1 to 9 such that in all rows, all columns and in nine 3X3 blocks every digit appears exactly once, while 
respecting a set of previously given digits in some of the cells (the so-called clues). Sudoku is an exact cover type 
constraint satisfaction problem^ and it is one of Karp's 21 NP-complete problems', when generalized to N X N 
grids'*. NP-complete problems are "intractable" (unless P = NP)^'' in the sense that all known algorithms that 
compute solutions to them do so in exponential worst-case time (in the number of variables N); in spite of the fact 
that if given a candidate solution, it takes only polynomial time to check its correctness. 

The intractability of NP-complete problems has important consequences, ranging from public-key cryp- 
tography to statistical mechanics. In the latter case, for the ground-state problem of Ising spin glasses (±1 spins), 
one needs to find the lowest energy configuration among all the 2" possible spin configurations, where N is the 
number of spins. Additionally, to describe the statistical behavior of such Ising spin models, one has to compute 
the partition function, which is a sum over all the 2" configurations. Barahona'', then IstraU' have shown that for 
non-planar crystalline lattices, the ground-state problem and computing the partition function are NP-complete'. 
Since there is little hope in providing polynomial time algorithms for NP-complete problems, the focus shifted 
towards understanding the nature of the complexity forbidding fast solutions to these problems. There has been 
considerable work in this direction, especially for the boolean satisfiability problem SAT (or fc-SAT), which is NP- 
complete for fc s 3. Completeness means that all problems in NP (hence Sudoku as well), can be translated in 
polynomial time and formulated as a A:- SAT problem, as shown for the first time by Cook and Levin^. Namely, any 
problem in NP can be solved via a small number of calls to a fc-SAT solver and a polynomial number of steps (in 
the size of the input) outside the subroutine invoking the fc-SAT solver. 

In fc-SAT we are given N boolean variables to which we need to assign Is or Os (TRUE or FALSE) such that a 
given set of clauses in conjunctive normal form, each containing k or fewer literals (literal: a boolean variable or its 
negation) are all satisfied, i.e., evaluate to TRUE. Just as for the spin glass model, here we also have exponentially 
many (2") configurations or assignments to search. 

In the following we treat algorithms as dynamical systems. An algorithm is a finite set of instructions acting in 
some state space, applied iteratively from an initial state until an end state is reached. For example, the simplest 
algorithm for the Ising model ground state problem, or the 3-SAT problem would be exhaustively testing 
potentially all the 2" configurations, which quickly becomes forbidding with increasing N. To improve perform- 
ance, algorithms have become more sophisticated by exploiting the structure of the problem (of the state space). 
Accordingly, now 3-SAT can be solved by a deterministic algorithm with an upper bound of 0(1.473*0 steps'*. 
Here we will only deal with deterministic algorithms that is, once an initial state is given, the "trajectory" of the 
dynamical system is uniquely determined. Thus, we expect that the dynamics of those algorithms that exploit the 
structure of hard problems will reflect the complexity inherent in the problem itself. Complex behavior by 
deterministic dynamical systems is coined chaos in the literature' and thus the behavior of algorithms for 
hard problems is expected to appear highly irregular or chaotic'"*. 
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Figure 1 | Sudoku and its boolean representation, (a) a typical puzzle with bold digits as clues (givens) . (b) Setup of the boolean representation in a 9 X 9 
X 9 grid, (c) Layer L4 of the puzzle (the one containing the digit 4) with 1-s in the location of the clues and the regions blocked out for digit 4 by the 
presence of the clues (shaded area). 



Although the theory of nonlinear dynamical systems and chaos is 
well-established, it has not yet been exploited in the context of optim- 
ization algorithms. One of the difficulties lies with the fact that most 
optimization algorithms are discrete and not easily cast in forms 
amenable to chaos theory methods. Recently, however, we have pro- 
vided'^ a deterministic continuous-time solver for the boolean satis- 
fiability problem fc-SAT using coupled ordinary differential equa- 
tions (ODE) with a one-to-one correspondence between the A:-SAT 
solution clusters and the attractors of the corresponding system of 
ODEs. This continuous-time dynamical system (CTDS) is in a form 
naturally suited for chaos theory methods, and thus it allows us to 
study the relationship between optimization hardness and chaotic 
behavior. Here we will focus only on solvable (satisfiable) instances, 
and thus the observed chaotic behavior will necessarily be tran- 
sient" We need to emphasize, however, that the dynamical 
properties characterize both the problem and the algorithm itself. For 
this reason, one compares the dynamical properties across problems 
of varying hardness using the same algorithm. Nevertheless, since 
there are problem instances that are hard for all known algorithms, 
the appearance of transient chaos with long lifetime should be a 
universal feature of hard problems. It is also important to observe 
that transient chaos is not an asymptotic behavior (where Nis 
the size of the problem in terms of input or number of variables), but 
it appears for finite N, and thus measures of chaos can be used to 
characterize and categorize the hardness of individual instances of 
finite problems. To illustrate this, here we first map the popular 9X9 
(hence finite) version of Sudoku into fc-SAT, then we solve it using 
our deterministic continuous-time solver'^. By analyzing the beha- 
vior of the corresponding trajectories of the CTDS we show the 
appearance of transient chaos when increasing the hardness of the 
Sudoku problems, and show that the level of hardness (taken from 
human ratings of the puzzles) correlates well with a chaotic invariant, 
namely the lifetime of chaos estimated as /c"', where ic is called the 
escape rate". We conclude with a discussion on algorithmic per- 
formance, dynamical properties and problem complexity. 

Results 

Sudoku as fc-SAT. Because our continuous-time dynamical system'^ 
was designed to solve fc-SAT formulae in conjunctive normal form 
(CNF), we first briefly describe how Sudoku can be interpreted as a 
+ l-in-9-SAT formula, and then how it is transformed into the 
standard CNF form. Further details are shown in the Methods 
section. 

In a Sudoku puzzle we are given a square grid with 9X9 = 81 cells, 
each to be filled with one of nine symbols (digits) D,je{l, ... ,9}, i,j 



= 1, . . ., 9 (with the upper-left corner of the puzzle corresponding to 
i = 1, j = 1). When the puzzle is completed, each of the columns, 
rows and 3X3 sub-grids (blocks partitioned by bold lines. Fig. la) 
must contain all the 9 symbols. Equivalently, all 9 symbols must 
appear once and only once in each row, column and 3X3 sub-grid. 

To formulate Sudoku as a constraint satisfaction problem (CSP) 
using boolean variables, we associate to each symbol (digit) an 
ordered set of 9 boolean variables (TRUE = "r', FALSE="0"). The 
digit Djj in cell (/, j) wiU be represented as the ordered set (xj, . . . ,xp 
with xJJe{0,1}, a = 1, 9, such that always one and only one of 
them is 1 (TRUE). In particular, = a is set to be equivalent to 
xfj = 5a^b, where Sa,b is the Kronecker delta function. This way we 
have in total 9 X 9 X 9 = 729 boolean variables xf^, which we can 
picture as being placed on a 3D grid (Fig. lb), with a corresponding 
to the grid index along the vertical direction, and hence a is the digit 
that is filling the corresponding {i, j) cell in the original puzzle. The 
corresponding horizontal 9X9 2D layer at height a will be denoted 
by La- Introducing the notion of such (horizontal) layers makes it 
easier to express the constraints of the Sudoku rules on its repres- 
entation by 0-s and 1 -s as described below. For example in the puzzle 
shown in Fig. la, we have Di ^ = 4. In the given vertical column the 
variable in the a* cell (that is in layer LJ is 9 = (5a 4 (shown as the 
boolean variable 1 filling the cell next to the digit 4 shown in red, in 
Fig lb). This setup allows us to encode the Sudoku constraints in a 
simple manner. They come from: 1) uniqueness of the digits in aU the 
(;', j) Sudoku cells, 2) a digit must occur once and only once in each 
row, column and in each of the nine 3X3 subgrids, and 3) obeying 
the clues. Constraint type 1) was already expressed above, namely 
that for every cell (;, j), in the set (xj, . . . ,x?) one and only one 
variable is TRUE, all others must be FALSE. Type 2) constraints 
are similar, e.g., in row / and layer La the set (x"j, . . . ,xfi^) must 
contain one and only one TRUE variable, all others must be false 
and this must hold for all rows and layers, etc. Observe that all 
constraints are in the form of a set of 9 boolean variables of which 
we demand that one and only one of them be TRUE, all others 
FALSE. When this is satisfied, we say that the constraint itself (or 
"clause") is satisfied, or TRUE. Such CSPs are called +l-in-fc-SAT 
and they are part of so-called "locked occupation problems", which is 
a class of exceptionally hard CSPs'"'''-". Type 3) constraints are gen- 
erated by the clues (or givens) which are symbols already filled in 
some of the cells and their number and positioning determines the 
difficulty of the puzzle. They are also set in a way to guarantee a 
unique solution to the whole puzzle. If there are given d clues, then 
this implies setting d boolean variables to TRUE, which means elim- 
inating exactly 4d constraints of type 1) and 2) (one vertical or 
uniqueness constraint, one row, one column and one 3X3 subgrid 
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constraint). Thus, Sudoku is a +l-in-9-SAT type CSP with N boo- 
lean variables and 324 — 4ci constraints of + l-in-fc- SAT type (A: £ 9). 
N is a complicated function of the positioning of the clues. The 
example in Fig. la has d = 22 clues with N = 232 unknown boolean 
variables. In layer L4 as illustrated in Fig. Ic there are 28 unknown 
boolean variables (white cells). These 28 variables appear in a total of 
17 constraints of + l-in-fc-SAT type. More precisely there is one + 1- 
in-2-SAT, six +l-in-4-SAT, four +l-in-5-SAT and sbc +l-in-6- 
SAT type of constraints related to L4. The other layers generate the 
remaining 324 — 4d — 17 = 219 of +l-in-/c-SAT type constraints 
(with k < 9). 

Since our continuous-time SAT solver has been designed to solve 
boolean satisfiability problems in conjunctive normal form (CNF), 
we need to bring the +l-in-fc-SAT type problems above into this 
form, and thus formulate it as a fc-SAT CNF problem. The CNF is a 
conjunction (AND, denoted by A) of clauses each clause expressed as 
the disjunction (OR, denoted by v) of literals. For fc-SAT in CNF 
there are N boolean variables Xj = {0, 1} and an instance is given as a 
propositional formula JF, which is the conjunction of M clauses Cm'- 
J-" = Ci A ■ ■ ■ A Cm A ■ ■ ■ A Ca(, with Cm = Zm, V . . . V z^^^ , ^ k 
and Zj is a literal, that is Zje^Xj, %j} {xj is the negation of x^). According 
to a well known theorem of propositional calculus, all boolean pro- 
positions can be formulated in CNF using De Morgan's laws and the 
distributive law, and thus any + l-in-/c-SAT type constraint as well. 
The Methods section describes how to translate + 1-in-fc-SAT type 
constraints into CNF. 

Once the transformation to CNF is completed we are left with N 
variables and M clauses of the type described above, called CNF clauses 
from here on. In our case the number of variables appearing in a CNF 
constraint has the property 1 < fc„ < 9. The parameters N, M and 
{fcm}m = i depend on the clues that are difficult to express analytically, 
but easy to determine computationally. The puzzle from Fig. la is 
ultimately formulated as a CNF SAT problem with N = 232 variables 
and M = 1718 CNF clauses. An often used parameter of a satisfiability 
problem is the number of CNF constraints per variable, or constraint 
density, y. = M/N, also used as a typical hardness indicator, however, as 
we show below, this is not an accurate measure of hardness. 

The continuous-time deterministic fc-SAT solver. In Ref [15] a 
continuous-time deterministic solver was introduced to solve fc-SAT 
problems in conjunctive normal form. The set of clauses specifying the 
constraints are translated into an M X N matrix: C = {c„„} with c„„- = 1 
if the variable Xj is present in clause m in direct (non-negated) form, 
namely XjeCm, c^i = — 1 if XjeCm and c^, = 0 if x, and Xi are both 
absent from C^. To every variable one associates a continuous spin 
variable s,e[ — 1,1] such that when s, = ±1 thenx; = (1 -1-5;) /2e{0,l}, 
and to every clause one associates the function: 

K„{s)=2-^"' n^{l-c„jSj), me{l,...,M}. (l) 

WehaveKmE[0,l] foraIlsE[ — 1,1]'^. It is easy to check that = Oonly 
for those Sj6{ — 1,-1-1} values for which the corresponding %,-s satisfy 
clause C„ (otherwise we always have > 0). That is, plays the role 
of an energy function for clause C„, and its ground state value of = 0 
is reached if C„ is TRUE, and only then. We also need the quantities , 
= -Km/(1 ~ CmiS,) that is, with the /-th term missing from the product in 
(1). Clearly, K,„,e[0,l/2]. The continuous time dynamical system 
introduced in'"" is defined via the set of {N + M) ordinary differential 
equations (ODEs): 

ds ** 

= ^ 2a„c„iK^i{s)K„,{s), i=l, . . . ,N (2) 

m— 1 

^=flmK,„(s), m = l,...M, (3) 



with the only requirements that s,(0)e[— 1,1], V; and a^{0) > 0, Vm. 
The latter implies from (3) that a^it) > 0, Vm, t. It was shown in Ref 
[15] that system (2-3) always finds the solutions to fc-SAT problems 
(encoded via the C matrix), when they exist, from almost all initial 
conditions (the exception being a set of Lebesgue measure zero). 
Here we give an intuitive picture for why that is the case. Due to (3) 
the auxiliary variables grow exponentially at rate K^. That is, the 
further is from its ground-state value of 0, the faster grows (in 
that instant). Moreover, the longer has been away from zero, the 
larger is a^, as seen from the formal solution to (3): a^{t) =0^(0) 

exp dzKmJ. Equation (2) can equivalently be written as a gradient 

descent on an energy landscape V(s, a), that is ds/dt = —V^V, where Vj 
is the gradient operator in the spin variables and V{s,a) = 
J2m '^mKl,{s). Clearly, y > 0 and V = 0 if and only if s is a 
fc-SAT solution, i.e., satisfies all the clauses (K^{s) = 0, Vm). 

From the behavior of the fl,„ variables discussed above it also 
follows that the least satisfied constraints will dominate V (terms 
with the largest flm'S). Without restricting generality, let the aiKf 
term be the most dominant at t. Then keeping only the dominant 
term on the rhs of (2) for those ; for which Cu 7^ 0 we get dsjdt = 
2fliCi,(l ^ Ci,s,)(_K'i,)^ or, equivalently: d(l — CuS^ldt = —(1 — Ci,s,) 
2fli (KijY. This shows that the term (1 — Ci,s,) is driven exponentially 
fast towards zero, that is towards satisfying (and all the other 
constraints containing this term). As Ki decreases, some other con- 
straint becomes dominant, and thus, in a continuous fashion, all 
constraints are driven towards satisfiability. The exponential growth 
guarantees that the trajectory is always pulled out of any potential 
well. When the problem is unsatisfiable, the system generates a cha- 
otic dynamics in [ — 1, 1]", indefinitely. For more details about the 
properties of the continuous-time dynamical system (CTDS) (2-3) 
see Ref [15]. 

Puzzle hardness as transient chaotic dynamics. Since Sudoku 
puzzles always have a solution, the corresponding boolean SAT 
CNF formulation also has a solution, and system (2-3) wUl always 
find it. The nature of the dynamics, however, will depend on the 
hardness of the puzzle as we describe next. 

In Fig.2a we show an easy puzzle with 34 clues (black numbers)''". 
After transforming this problem into SAT CNF, we obtain N = 126 
and M = 717, with a constraint density of a = M/N = 5.69. As 
described above, in our implementation there is a continuous spin 
variable associated to every boolean variable x"- in every 3D cell {i,j, 
a). In the right panels of Fig. 2 we show the dynamics of the spin 
variables in the cells of the 3X3 grid formed by rows 4-6 and 
columns 7-9. The sj(t) curves are colored by the digit a they repres- 
ent (a = 1, 9) as indicated in the color legend of Fig 2. The 
dynamics was started from a random initial condition. Indeed, our 
solver finds the solution very quickly, in about 15 time units, for the 
easy puzzle in Fig.2a. 

In Fig. 2b we show the dynamical evolution of variables for an 
ultra-hard Sudoku instance with only 21 clues. This puzzle has been 
listed as one of the world's hardest Sudokus, and even has a special 
name: "Platinum Blonde"^'"^^, and it was the most "difficult" for our 
solver among all the puzzles we tried. After transforming it into SAT 
CNF, we obtain N = 257 variables and M = 2085 constraints. Not 
only that we have twice as many unknown variables but the con- 
straint density a = M/N = 8.11 is also larger than in the previous 
case, signaling the hardness of the corresponding SAT instance. The 
complexity of the dynamics in this case is seen in the right panel of 
Fig. 2b, exhibiting long chaotic transients before the solution is found 
at around t ~ 150. For an animation of the dynamics for a similarly 
hard puzzle" see Ref [24] . 

We can also observe from the right panels in Fig 2 that there is one 
dominating digit (a-value), corresponding to which vertical cell at 
that given {i,j) grid cell has the largest s"j value. This can be taken as 
the digit Dy the solver is considering in the given grid cell {i,j) at that 
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Figure 2 | Solving Sudoku puzzles with the deterministic continuous-time solver (2-3). (a) presents an easy puzzle with the evolution of the 
continuous-time dynamics shown within a 3 X 3 grid (rows 4-6, columns 7-9). (b) shows the same, but for a known, ultra-hard puzzle called "Platinum 
Blonde"^'"^'. The trajectories in the right panels show the evolution of the analog variables ( f ) colored by the corresponding digit a. Thus for each cell ( !, 
j) we have 9 such running trajectories, but they cannot always be discerned as many of them are running on top of each other, close to -1, as seen in (a). 



moment. We wiU use this observation to provide below an alternate 
illustration of the dynamics' transiently chaotic behavior. Let us fix a 
random initial condition except for two chosen variables that are 
varied along the points of a square grid within the domain [ — 1, 1]^. 
There is no particular relevance as to which pairs of variables are 
chosen to be varied, let us denote them by Si and Sj. Let us choose 
an arbitrary empty ceU (p, q) in the original Sudoku puzzle and mon- 
itor the dominating digit in it at time t. We wiU color the initial 
conditions in the plane (si, S2) according to the dominating digit in 
(p, q) at time t. This will provide a map expressing the "sensitivity to 
initial conditions" that varies across time. Since aU puzzles have solu- 
tions, the maps eventually assume one solid color according to the 
digit of the solution in the monitored cell, however, for hard puzzles, it 
may assume highly complex patterns before it does that, as shown in 
Fig. 3. Here we show these colormaps for the easy and hard Sudoku 
puzzles from Fig. 2 at times t = 10, 15, 20. For the easy puzzle (top 
row) the cell was chosen to be (p, g) = (1, 1). At time f = 10 the whole 
map shows Di i = 6 (orange), which is not the solution digit (it is still 
searching for the solution). At time t= 15, however, we see two clearly 
separated domains, in one of them Di 1 = 6, in the other Di 1 = 9 (light 
blue), the latter being the correct digit. As time passes, the orange 
(incorrect) domain shrinks, because trajectories from an increasing 
number of initial conditions find the solution. At t = 20 almost the 
whole map shows the correct digit = 9, except for a thin line. 

In the case of the hard Sudoku puzzle (bottom row in Fig. 3, (p, q) = 
(6, 8)) more colors enter the picture over time, in a fractal-like pattern. 
Eventually the whole map becomes dark blue (digit 4) corresponding 
to the solution in this cell. Changing the initial condition slightly about 
this fractal set (which is really a time-dependent fractal boundary) 
may result in a completely different digit (color) being considered in 
cell (p, q) at time t. This sensitivity to initial conditions is indicative of 
the chaotic behavior of the (deterministic) search dynamics. 

The appearance of transient chaos is a fundamental feature of 
the search dynamics and can be used to separate problems by their 



hardness. In Ref [15] we have shown that within the thermodynamic 
limit {N^ 00, 00, a = M/N = const.) of random fc-SAT ensem- 
bles this appears as a phase transition at the so-caUed chaotic trans- 
ition point V-y in terms of the constraint density a = M/N. Since there 
is no "thermodynamic limit" for 9X9 Sudoku problems {N < 729), 
one cannot define a simple order-parameter and use it to rate prob- 
lem hardness in the same way'^. However, once a problem is given, 
the corresponding dynamical system (2-3) is well defined, and so is 
its dynamical behavior. Even though we do not have a weU-defined 
ensemble-based statistical order parameter (which has little meaning 
for specific SAT instances anyway), here we show next how can we 
use a well-known invariant quantity from non-linear dynamical sys- 
tem's theory to categorize problem hardness for specific instances. 

A Richter-type scale for Sudoku hardness. As suggested by the two 
examples in Fig. 2, the hardness of Sudoku puzzles correlates with the 
length of chaotic transients. A consistent way to characterize these 
chaotic transients is to plot the distribution of their lifetime. Starting 
trajectories from many random initial conditions, letp(t) indicate the 
probability that the dynamics has not found the solution by time t. A 
characteristic property of transient chaos"''^ in hyperbolic 
dynamical systems is that p{t) shows an asymptotic exponential 
decay: p{t) ~ e "', where k is called the escape rate. For an easy to 
read text on transient chaos see the book". The escape rate, a 
measurable quantity, theoretically can be expressed as a zero of the 
spectral determinant of the evolution operator corresponding to the 
dynamical system (2-3) and well approximated using the machinery 
of cycle expansions based on dynamical zeta functions'^. It is an 
invariant measure of the dynamics in the sense that it characterizes 
solely the chaotic nonattracting set in the phase space of the system, 
and it does not depend on the distribution of the initial conditions, its 
support, or the details of the region from where the escape is 
measured (as long as it contains the non-attracting set)". 

In Fig. 4a we plot the distribution p{t) in log-linear scale for several 
puzzles gathered from the literature. The distributions were obtained 



SCIENTIFIC REPORTS | 2 : 725 | DOI: 1 0. 1 038/srep00725 



4 



t = 10 t = 15 t = 20 




Figure 3 | Puzzle hardness as chaotic dynamics. We color the points of a 10' X 10' grid in an arbitraryplane (si, S2) at time instant taccordingto the digit 
Dp^ the solver is considering in an arbitrary but fixed cell (p, q) at that instant, given that we started the trajectory of the CTDS from those grid-points. For 
these initial conditions only the points in the (si, S2) plane were varied, aU other spin values were kept fixed at the same randomly chosen values. For the 
same easy problem as in Fig. 2, and for (p, q) = ( 1 , 1 ) (top row of panels) almost aU initial conditions in this plane involve only two digits, and after t = 20 
the corresponding trajectories have converged to the solution digit (9, light blue), except for a thin line, which, however, will also become light blue. 
The bottom row of panels shows the same for the hard problem of Fig. 2 based on what happens in the cell (p, q) = (6, 8). The strong sensitivity to initial 
conditions appears as fractal structures of increasing complexity as time goes on, before eventually everything converges to the same color/digit (dark 
blue, corresponding to digit 4, not shown). 



from over 10* random initial conditions. The decay shows a wide 
range of variation between the puzzles. For easy puzzles the transi- 
ents are very short, p{t) decays fast resulting in large escape rates but 
for hard puzzles k can be very small. Fig. 4b shows a zoom onto the 
p{t) of hard puzzles. For example, for the puzzle in Fig.2a we obtain /c 
= 0.156, whereas for Fig.2b (Platinum Blonde) the escape rate is k = 
0.00026. In spite of the large variability of the decay rates, we see that 
in aU cases the escape is exponentially fast (the faster than exponen- 
tial appearance for some of the curves in Figures 4a,b is due to finite 
size statistics and the finiteness of the time interval considered). 

The several orders of magnitude variability of k naturally suggests 
the use of a logarithmic measure of k for puzzle hardness, see Fig.4c, 
which shows the escape rates on a semilog scale as function of the 
number of clues, d. Thus, the escape rate can be used to define a kind 
of "Richter"-type scale for Sudoku hardness: 

'/= -logio('^) (4) 

with easy puzzles falling in the range 0 < £ 1, medium ones in 1 < 
»7 s 2, hard ones in 2 < £ 3 and for ultra-hard puzzles rj > 3. 

We chose several instances from the "Sudoku of the Day" web- 
site^" in four of the categories defined there: easy (black square), 
medium (red circle), hard (green X) and absurd (blue star). These 
ratings on the website try to estimate the hardness of puzzles when 
solved by humans. These ratings correlate very well with our hard- 
ness measure )}, giving an average hardness value of (rj) = 0.816 for 
easy, {t]} = 1.439 medium, {r]) = 1.782 for hard and (;]) = 1.809 for 
what they call absurd. Another site we analyzed puzzles from is 
"Extreme Sudoku"^^ (brown + signs on Fig.4). It claims to offer 
extremely hard Sudoku puzzles, their categories being: evil, excessive, 
egregious, excruciating and extreme. Indeed those puzzles are 



difficult with a range of j; e [1.1, 1.9] on the hardness scale, however, 
still far from the hardest puzzles we have found in the literature. 
Occasionally, daily newspapers present puzzles claimed to be the 
hardest Sudoku puzzles of the year. In particular, the escape rate 
for the Caveman Circus 2009 winner^' (turquoise diamond) and 
the Guardian 2010 hardest puzzle^^ (maroon diamond) are indeed 
one order of magnitude smaller than the hardest puzzles on the daily 
Sudoku websites, placing them at i; = 2.93 and )} = 2.82 on the 
hardness scale. The USA Today 2006 hardest puzzle^", however, does 
not seem to be that hard for our algorithm having = 2.17 (magenta 
diamond). Eppstein^' gives two Sudoku examples (orange left-point- 
ing triangles) while describing his algorithm, one with i] = 1.288 and 
a much harder one with )} = 2.017. Elser et aV* present an extremely 
hard Sudoku (black filled circle), which has an escape rate of k = 
0.0023 resulting in 17 = 2.639. 

The smallest escape rates we have found are for the Sudokus listed 
as the hardest on Wikipedia^^'^'' (red triangles). The five puzzles, 
which we tested are called Platinum Blonde, Golden Nugget, Red 
Dwarf, coly013 and tarx0134. They have a hardness in the range 3 < 
f; < 3.6, the ultra-hard Platinum Blonde (shown in Fig.2b) being the 
hardest with 1? = 3.5789 (k = 0.00026). 

While the escape rate correlates surprisingly well with human 
ratings of Sudoku hardness, it is natural to expect a correlation with 
the number of clues, d. Indeed, as a general rule of thumb, the fewer 
clues are given, the harder the puzzle, however, this is not universally 
true'. Here we tested a few instances with minimaP", that is 17 clues 
and almost minimal 18 clues (orange filled circles)'' As seen from 
Fig.4c, these are actually easier (1.2 < i] < 2.4) than the hardest 
instances with more d = 21, 22 clues. In Fig.4d we then plot the 
escape rate as function of the constraint density a = M/N, leading to 
practically the same conclusion. This is because the constraint 
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Figure 4 | Escape rate as hardness indicator, (a) shows the distribution in log-linear scale of the fraction p( f) of 10* randomly started trajectories of (2-3) 
that have not yet found a solution by analog time tfor a number of Sudoku puzzles taken from the literature (see legend and text) with a wide range of 
human difficulty ratings. The escape rate is obtained from the best fit to the tail of the distributions, (b) is a magnification of (a) for hard puzzles, (c) and 
(d) show the escape rate k in semilog scale vs the number of clues d and constraint density a indicating good correlations with human ratings (color 
bands), (e) shows the relationship between the number of clues d and a for the puzzles considered. 



density a is essentially linearly correlated with the number of givens 
d, as shown in Fig.4e. The apparent non-monotonic behavior of 
puzzle hardness with the number of givens or constraint density is 
due to the fact that hardness cannot simply be characterized by a 
global, static variable such as d or a, but it also depends on the 
positioning pattern of the clues'. 

Discussion 

Using the world of Sudoku puzzles, here we have presented further 
evidence that optimization hardness translates into complex dynam- 
ical behavior by an algorithm searching for solutions in an optimal 
fashion. Namely, there seems to be a trade-off between algorithmic 
performance and the complexity of the algorithm and/or its beha- 
vior. Simple, sequential search algorithms have a trivial description 
and simple dynamics, but an abysmal worst-case performance (2''), 
whereas algorithms that are among the best performers are complex 
in their description (instruction-list) and/or behavior (dynamics). 
This happens because in order to improve performance, algorithms 



have to exploit the structure of the problem one way or another. As 
hard problems have complex structures, the dynamics of the algo- 
rithms should be indicative of the problem's hardness. However, as a 
word of caution, observing complex dynamics performed by some 
black-box algorithm does not necessarily imply problem hardness. 
For example, one could consider any dynamical system that is guar- 
anteed to visit all the 2" discrete states/configurations but perhaps 
with some arbitrary, complex behavior. This algorithm will always 
find solutions if they exist (by checking for satisfiability after every 
new state visited). But its instruction list would have no relevance to 
the problem itself (apart from the checking instructions to see if the 
new state satisfies the problem) and thus, it could take long times to 
find solutions even for problems that are otherwise easily solved by 
other algorithms. Hence, dynamical properties can only be regarded 
as descriptors of problem hardness if they are generated by algo- 
rithms that: 1) exploit the structure of the state space of the problem 
and 2) they show similar or better performance compared to other 
algorithms on the same problems. 
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The continuous-time dynamical system'^ (2-3) as a deterministic 
algorithm does have these features: 1) the search happens on an 
energy landscape V= 5Zm '''"^m that incorporates simultaneously 
all the constraints (problem structure) 2) it solves easy problems 
efficiently (polynomial time, both analog and discrete) and 3) it 
guarantees to find solutions to hard problems even for solvable cases 
where many other algorithms fail. Although it is not a polynomial 
cost algorithm, it seems to find solutions in continuous-time t that 
scales polynomially with N'"'. These features and the fact that the 
algorithm is formulated as a deterministic dynamical system with 
continuous variables, allows us to apply the theory of nonlinear 
dynamical systems on CTDS (2-3) to characterize the hardness of 
boolean satisfiability problems. In particular, via the measurable 
escape rate k, or its negative log-value we can provide a single- 
scalar measure of hardness, well defined for my finite instance. We 
have illustrated this here on Sudoku puzzles, but the analysis can be 
repeated on any other ensemble from NP. Having a mathematically 
well-defined number to characterize optimization hardness for spe- 
cific problems in NP provides more information than the poly- 
nomial/exponential-time solvability classification, or knowing what 
the constraint density a = MIN is (the latter being a non-dynamic/ 
static measure). Moreover, within the framework of CTDS (2-3), 
dynamical systems and chaos theory methods can now be brou- 
ght forth to help develop a novel understanding of optimization 
hardness. 

Methods 

Type 1) constraints (main text) impose the uniqueness of the symbol Dy in a given 
cell, expressed as a -t- l-in-9-SAT constraint on the vector: 

(4,4... 4). (5) 

that is, one and only one component in the boolean vector above can be TRUE (or 1), 
all others must be FALSE {or 0). Having 9X9 cells in the puzzle, this gives in total 81, 
+ l-in-9-SAT constraints. 

Type 2) constraints on rows, columns and sub-grids further impose that in every 
layer we have the following 27, +l-in-9-SAT constraints: 



clause defined on the (71, 72* -■■>7fc) variables can be written as one fc-SAT andA:(A;-l)/2 
of 2-SAT constraints: 



Rows : 

(4,4,..., 4), i=l,...,9 

Columns : 

(x«.,4,...,4)j=l,...,9 



(6) 



(7) 



Subgrids : 



■^m 2,/i -M '-^m -I- 2,/j 2 ■'^m -I- 2,« -I- 3 ' 
^ 



(8) 



— 0,3,( 



with the same meaning as before, that is in all the boolean vectors above one and only 
one variable can be TRUE all others must be FALSE. Together with the 8 1 constraints 
of type 1) we thus have a total of 9 X 27 + 81 — 324 constraints in + l-in-9-SATform. 

Finally, type 3) constraints are imposed via d given digits or clues. It was only 
recently shown that uniqueness of a solution demands that d ^ 17^". As discussed in 
the main text, each clue will eliminate 4 constraints: in its vertical tower, its column, its 
row and the 3X3 sub-grid containing the clue. For example, let us examine layer L4 
(Fig.lc) of the puzzle shown in Fig. la. There are three clues of4 in cells (1,9), (3, 3), (4, 
4) and thus ^ — I, Xj ^ — I, x"^ ^ — I have to be fixed as TRUE in I4. In order to satisfy 
the constraints, the other variables in the same rows, columns, blocks and vertical 
columns must be set to FALSE. The unknown variables left in the SAT problem will be 
those in the light cells of Fig.lc. (The other clues will eliminate constraints and 
variables in other layers and vertical columns.) The total number of unknown vari- 
ables N depends on d and on the placement of clues. The number of constraints is 
always 324 - 4d, however the number of variables in a clause can vary. For example in 
Fig. 1 c the constraint corresponding to the second row in layer L4 has only 2 unknown 
variables left (+ l-in-2-SAT). 

After the unknown boolean variables and the constraints have been identified we 
need to transform the formula into CNF. This is done easily, as any + l-in-fc-SAT 



(yi Aiyi'^yj) 



(9) 



The disjunction (V) of the first k variables enforces that at least one variable must be 
true, but the rest of {k{k - 1)/2) 2-SAT type constraints ensure that only one of them is 
allowed to be true. 
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